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We investigate the miscible Rayleigh- Taylor (RT) instability in both 2 and 3 dimensions 
using direct numerical simulations, where the working fluid is assumed incompressible 
under the Boussinesq approximation. We first consider the case of randomly perturbed 
interfaces. With a variety of diagnostics, we develop a physical picture for the detailed 
temporal development of the mixed layer: We identify three distinct evolutionary phases 
in the development of the mixed layer, which can be related to detailed variations in 
the growth of the mixing zone. Our analysis provides an explanation for the observed 
differences between two and three-dimensional RT instability; the analysis also leads us to 
concentrate on the RT models which (1) work equally well for both laminar and turbulent 
flows, and (2) do not depend on turbulent scaling within the mixing layer between fluids. 
These candidate RT models are based on point sources within bubbles (or plumes) and 
interaction with each other (or the background flow) . With this motivation, we examine 
the evolution of single plumes, and relate our numerical results (of single plumes) to a 
simple analytical model for plume evolution. 



1. Introduction 

The phenomenology of mixing due to the Rayleigh- Taylor instability, the instability 
of an in terface separating fluids of different densities subject to gravity ( |Chandrasekhar 



1961)1 and references therein), can be summarized as follows: Bubbles (spikes) of lighter 
(heavier) fluid penetrate into the heavier (lighter) fluid, and leave behind them a region 
of mixed fluid. As the instability enters the nonlinear regime, fluid motions in this mixing 
region become highly irregular. The envelope, or edge, of the mixing region is observed 
to be deflned by fast large-scale motion, which tend to be dominated by the merging of 
expanding smaller bubbles (spikes). The irregular fluid motions within the mixing zone 
are often regarded as "turbulent" in the sense that the flow is chaotic in the wake of the 
bubbles (spikes). 

A principal focus of experimental and theoretical study has been the general prop- 
erties of this mixing zone, whose broadening in time is commonly characterized by the 
"envelope velocity" (Ve) and the penetration lengths /it, and hg (for bubbles and spikes, 
respectively). Using simple analytical models for the interpenetration of the two flu- 



ids, one can show (Youngs (1984) and references therein) that in the nonlinear regime 
Ve is proportional to time (t), and that the penetration depths hb^s are proportional 
to and depend linearly on the gravitational acceleration g and the Atwood number A 
[= (pi — p2)/{pi +P2), where pi,2 is the density of the heavier (lighter) fluid, respectively], 
i.e.. 
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l'b,s 



(1.1) 



a can be thought of as a measure of the efficiency of potential energy release; experimental 
mea surements of a give results in th e range of 0.03 ^ 0.06 (s ee Sharp (1984) , Read (1984) 
and Youngs (1984) ). Most recently, Schneider et al. (1998) show that a lies in the range 
of 0.05 ~ 0.06 for an Atwood number of 0.34; as the Reynolds number is high 10^) 
in these Linear Electric Motor (LEM) experiments, the instability enters the nonlinear 
regime within several e-folding times {'^ 0.9 ms), so that the above scaling law seems to 
hold from the time that the first few measurements of mixing zone width are made. 

There are many possible reasons for the observed variation of a found in the literature, 
when scaling of the form Eq. (^]^) is observed during the experiment or simulation. Set- 
ting aside problems such as incompletely controlled experiments or insufficiently resolved 
computations, it is important to establish whether the nonlinear (long-time) evolution 
of the instability is sensitive to details in the experiment such initial conditions; and to 
what extent one can really talk about universal scaling of the mixed layer widt h. These 
questions have been addressed to some extent in the literature. For example. Youngs 



(1991)1 has studied the variation of a with time during the course of R-T evolution; he 
reported that larger values of a ~ 0.06 occur during the early phase of the 3-D prob- 
lem, before any development of small-scale turbulence; however, as small-scale motions 
develop, Youngs (1991) reports that the 3-D growth rate slows to that characteristic of 



the 2-D case [a ^ 0.04). It is also found that a may depend on the Atwood number A as 
A ^ 1 (G. Dimonte, private communication). In a variation on this theme, it has been 
observed in simulations that a can depend on dimensionality: a for three-dimensional cal- 



culations is found to be larger than from 2-D simulations. Sakagami & Nishihara (1990) 



have found a for three-dimensional calculations to be at least 4 times higher than for 2-D 
simulations (in spherical systems); and more detailed phenomena can be found in Yab^ 
(1991)1 , |Li (1993)1 and |Hc (1999)| . [Town fc Bell (1991)| found larger values for a in 3-D only 
during the early nonlinear stage, and showed that as small scale 3-D turbulence devel- 



oped, the 3-D growth slows to that of the 2-D case, similar to results reported by Youngs 



(1991). Other 3-D numerical simulations, including investigations of single-mode initial 
perturbation (Dahlburg & Gardner (1990), Tryggvason& Unverdi (1990), Kane & Arnett 
(1998)) and more general initial perturbations (Cook (1998)), obtain a range of values 
for a generally restricted to 0.05 ~ 0.06. In this paper we also explore the dependence 
of the RT instability on dimensionality. We focus on quantitative comparisons between 
two and three dimensions, and these comparisons provide more complete insights to the 
understanding of RT instability. In the case of experiments conducted with gases, it has 
also been shown that the interface perturba tion growth can b e measurably affected by 
mass diffusion effects between the two layers Duff et al. (1962) . These experimental and 
computational results argue strongly that the precise value of a is sensitive to a variety 
of details in the experiments, ranging from the specific nature of the initial state to the 
dimensionality of the perturbation. It is the aim of this paper to initiate a study of this 
issue: In this paper we consider some of the effects which influence the actual value of 
the scaling coefficient a as the Rayleigh- Taylor instability develops in time; a central 
question we address is whether a single scaling regime is always to be expected during 
the course of non- linear Rayleigh- Taylor development. 

In order to make our study more tractable, we have focused only on the low Atwood 
number regime. The advantage of this regime is that numerical techniques developed 
for incompressible flows can be readily applied to flows satisfying the Boussinesq equa- 
tions; such computations are substantially more economical than calculations of the fully 
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compressible problem, and thus allow a more extensive exploration of the parameter 
space other than Atwood number. Physically, our calculations correspond to, for exam- 
ple, cases in which the density contrast is introduced by the distribution of a scalar field 
(e.g., salt or temperature), for which the Atwood number is usually small because the 
density varies weakly with the scalar "concentration" . Species or temperature diffusion 
occurs in such systems, and is therefore included in our simulations; we are therefore 
considering the miscible version of the Rayleigh- Taylor instability; in this Boussinesq 
approximation, the Atwood number A and yet gA is a finite constant. An essential 
element in computations in this regime is that the calculations remain fully resolved 
throughout the evolution of the instability. This important constraint limits the range 
of allowable Prandtl number Pr (ratio of viscosity to thermal diffusivity) or Schmidt 
number Sc (ratio of viscosity to species diffusivity) i n ou r calculations; we typically take 



these ratios to be of order unity. [Tryggvason (1988)| and |Aref fc Tryggvason (1989)| have 



adopted similar ideas in their earlier studies of immiscible Rayleigh- Taylor instability; 
however, in their case, the density is not coupled to a scalar field, and the buoyancy term 
in the Boussinesq Navier-Stokes equations (induced by the sharp density gradient across 
the interface) is the only term related to the weak stratification. 



There are unfortunately few experimental studies of this regime: Linden et al. (1994) 
have studied the miscible Rayleigh- Taylor instability experimentally, and via simulations. 
In their experiments, they place a layer of brine above a layer of fresh water, separated 
by an aluminum barrier. The experiment is initiated by sliding the barrier horizontally 
through one sidewall; the Atwood numbers in their experiments are in the range of 10"'^ ~ 
10^^ (as the density contrast due to the brine concentration is small when compared to 
the ambient density) and the Schmidt number S'c is ^ 10'^. Their experimental estimates 



for a are around 0.044. Linden et al. (1994)| also performed numerical simulations for an 



Atwood number A — 0.2, and (numerical) Schmidt number of order unity. After imposing 
various perturbations at the interface, and determining a as a function of the ratio of the 
initial displacement of the fluid to the perturbation wavelength, they concluded that the 
smaller this ratio is, the higher the value of a. These authors argued that the removal 
of the barrier in the experiment corresponds to a long wavelength perturbation at the 
interface, and therefore the comparison between numerical simulation and experiments 
makes sense only if the long wavelength perturbation is included as part of the initial 
perturbation in the simulation. By implication, they thus suggest that the value of a 
obtained even in the long-time limit can depend on the initial conditions; in other words, 
these systems have long-term memory that is not erased by turbulent motions within the 
mixing interface. The dependence of RT instability on the initial conditions has recently 



been studied by Cook (2000) . In this paper, we focus on a study of the physical effects 
which determin a, and in particular investigate effects which lead to a departure from 
simple t^ scaling of the mixing zone depth. A detailed study of the scaling law itself is 
now in progress, and led by G. Dimonte (the a group). 

Our paper is organized as follows: We first formulate the problem, outline the numer- 
ics, and present evidence for the convergence and accuracy of our calculations. In the 
following section (§3) we summarize our results for multi-mode interface perturbations. 
We next discuss existing (analytical) models relevant to our studies (§ 4), which motivate 
a synthesis of a new model; the following section (§ 5) we explore these ideas in more 
detail via numerical studies of isolated plume evolution in both 2-D and 3-D. Our results 
are discussed in § 6, while our conclusions are presented in § 7. 



2. Formulation and Methods 
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2.1. Formulation of the problem 
We consider two vertically stacked fluid layers of different density, governed by the Boussi- 



nesq equations in a rectangular box (Cliandrasekhar (1961)). We adopt periodic bound 



ary conditions along the horizontal directions; the top and bottom boundaries are no- flux 
and no-slip. Under the Boussinesq approximation, density variations only appear in the 
buoyancy term, and are small when compared to the mean density. There are a vari- 
ety of ways of achieving a small density jump across the horizontal inerface: one may 
choose gasses of slightly different mean molecular weight, or a single fluid into which a 
'contaminant' is added (such as temperature, sugar or salt in water) which changes the 
fluid's density slightly, such vertical density variation depends linearly on temperature; 
colder (heavier) fluid sinks and warmer (lighter) fluid floats as gravity points downwards. 
Thus, in our formulation, the density jump across the interface is replaced by (for ex- 
ample) a temperature difference: colder fluid on top of warmer fluid. The equations for 
the miscible, incompressible, Rayleigh- Taylor instability in the Boussinesq limit are then 
essentially the incompressible Navier-Stokes equations coupled to an advection-diffusion 
equation for the 'contaminant', here temperature: 

dtu+{u-V)u = -VP + vV^u + PgTg, (2.1) 
V-M = 0, (2.2) 
dtT +{u-V)T = kV^T, (2.3) 

{v is the kinematic viscosity, (3 is the volumetric expansion coefficient, and k is the thermal 
diffusivity.) We remark that the initial, unperturbed temperature profile diffuses on the 
diffusion time scale. We further remark that one could alternatively treat the case of 
weak stratification via a vertically-imposed salinity (or any other type of concentration) 
gradient; in this way, one is able to explore a wide range of Prandtl/Schmidt numbers. 
There can be a variety of top and bottom boundary conditions for the scalars, such as 
fix values or fix fluxes at the boundaries. In our simulations, we adopt no-flux boundary 
conditions, which are typical in a run-down experiment. It is important to note here that 
since our equations allow for temperature diffusion, they are equivalent to systems in 
which the 'contaminant' is salt and allow for saline diffusion, or to systems describing 
two miscible fluids of slightly different density and allow for mass diffusion. 

In our simulations, we are interested only in regimes where the thermal diffusion time 
is much longer than the dynamical time scales of interest. The reason for this constraint 
is as follows: Unlike the case of immiscible Rayleigh- Taylor instability, the presence of 
thermal diffusion formally does not allow for a static background temperature profile prior 
to any perturbation. Nevertheless, as we perturb the temperature at the interface, we 
find that if the diffusivity is sufficiently small, then the initial growth of the perturbation 
is exponential and thus mimics the linear regime in standard stability analyses for which 
an initial static equilibrium exists. Thus, long thermal diffusion times are essential if 
the dynamics is to be dominated by the Rayleigh- Taylor instability. The same kind of 
argument applies in the saline case, or in the case where mass diffusion can occur. Thus, 



we do not expect to see diffusion effects as reported by Duff et al. (1962); and our results 
indeed show that at least in the linear regime, our growth rates are the same as in the 
inviscid nondiffusive case. 

Our problem is defined by a number of distinct spatial and temporal scales: the pertur- 
bation wavelength(s) A (for single-mode perturbation) or spatial perturbation spectrum 
(for multi-mode perturbation), the amplitude of the initial displacement of the fluid at 
the interface aA, and the initial interface thickness d. Two timescales are of most inter- 
ests here: the dynamical time defined by the free-fall time and the scalar diffusion time. 
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In the present study, we focus on cases where the dynamical time scale is much smaller 
than the diffusion time scale. In this limit, the diffusive interface between the hot and 
cold material remains more or less fixed at width d throughout the evolution; the com- 
putation's spatial resolution is then required to at least resolve this interface throughout 
the calculation. Furthermore, we impose perturbations of amplitude a such that a > d in 
order to capture the initial perturbation properly. Finally, the range of Atwood numbers 
we explore is consistent with values obtained in laboratory experiments for the miscible 
RT instability; and the ratio of viscosity to diffusivity is kept to values greater than or 
equal to 1. 

2.2. Numerics: methods and validation 

We have used two distinct numerical codes to solve the above equations, based on our 
desire to insure that our results are independent of the computational method used to 
solve equations (2.1)-(p!^). We have compared both the mixing zone extent and the 
evolution of fluid structures (obtained from these two codes) , as to determine the fldelity 
of our calculations: The two codes lead to results within expected range of errors for the 
same initial conditions and boundary conditions. 




Figure 1. Panel (a): Comparison between the spectral-element code (solid line) and the spectral 
code (dashed line) for simulation of single mode perturbation at the interface. Plotted is the 
mixing zone width /i as a function of time r. Panel (b): The spectral code shows exponential 
convergence for a fixed interface thickness of 0.15 in the single mode simulation. The error is 
defined as the deviation from the simulation with a resolution of 1024. 

The first code is based on the use of pseudospectral methods to simulate the RT in- 
stability. The spatial discretization is Fourier in the periodic, horizontal directions and 
Chebyshev in the gravitational direction. The temporal discretization is 3rd order Adams- 
Bashforth for the nonlinear terms and Crank-Nicholson for t he Laplacian t erms. We use 
Werne's tau-correction scheme to achieve incom pressibility ( Werne (1995) ); the code is 
parallelized using MPI (Dubey & Clune (1999)), and has been run on the Cray-T3E su- 
percomputer at PSC. The second code is based on the spectral element method ( [Tufo &: 
Fischer (1999)[ ). It combines the spectral accuracy /efficiency and the geometric flcxibil 



ity of finite elements, and has been widely applied to a variety of problems. Figure ^ (a) 
demonstrates the satisfactory agreement between the two codes. We plot the penetration 
depth h of the single mode perturbation as a function of time r. The solid line is from 
spectral-element code, and the dashed line is from spectral code. The spectral resolution 
is 256 by 512 and the spectral element code has an equival number of total grid points 
for this comparison. The average difference in h is as small as 10"''. By comparing results 
between spectral-element and spectral code, we have validated the usage of spectral code 
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N 


E{ti) 


errori 


error g 


7 


1.11498657 


0.003963 


0.313602 


9 


1.11519192 


0.003758 


0.001820 


11 


1.11910382 


0.000153 


0.004407 


13 


1.11896714 


0.000016 


0.000097 


15 


1.11895646 


0.000006 


0.000041 



Table 1. Spatial convergence of the spectral-element code for the Orr-Sommerfeld problem: 
K = 15, dt = .003125. E{t\) is the energy of the perturbation at ti = 25.1437, the one period of 
osciUation for the TS wave, errori is the error in energy when compared to theory, and errorg 
is the error in growth rate. 



in this particular situation, where the smaU transitional region from one fluid density 
to another may undermine the spectral convergence. In addition to the code-code com- 
parison, we have also performed standard convergence tests on the spectral code when 
applied to the miscible RT simulation to assure ourselves of the spectral convergence with 
the presence of a thin interface. We obtain the expected spectral convergence if the inter- 
face is smooth and well resolved. The ideal density discontinuity (or, in our case, a sharp 
temperature jump) cannot be resolved in spectral codes without special treatment; we 
therefore use a hyperbolic tangent or an error function for the background temperature 
field to make sure that the initial temperature field is smooth. From the smoothness of 
the interface, we obtain spectral convergence as long as we resolve the interface. Figure 
|l| (b) shows the spatial convergence of the spectral code when we fix the interface thick- 
ness and use single mode perturbation as the initial condition. The convergence test for 
the spectral element code is conducted separately for the Orr-Sommerfeld equation as 
shown in table ^ The convergence is non-monotonic due to the fact that the growth rates 
oscillate about the analytical value. However, spectral convergence is clearly attained. 

Having demonstrated the validation of the spectral code, we present results from sim- 
ulating the system with the spectral code in the rest of the paper. We adjust both the 
resolution and the interface thickness so that there are enough grid points across the 
background temperature profile. 

3. "Random" perturbations at the interface: 2-D versus 3-D 




L 

(*} I?) 

Figure 2. Initial random perturbation at the interface {z — 0.06). Panel (a) is the top view of 
the perturbation and panel (b) is the power spectrum of the perturbation. 
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In this section, we focus on the consequence of perturbing the interface with a spa- 
tiaUy randomized disturbance. The initial perturbation imposed at the interface in the 
3-D simulation is shown in figure |^. This perturbation is the outcome of phase random- 
ization of a spatial spectral distribution which peaks at fc ^ 8 and whose bandwidth is 
approximately 4 centered at the peak. This initial perturbation has been defined as one 
of the canonical test problems of the "alpha group" collaboration (A. Cook & G. Di- 
monte, private communication). The perturbation spectrum (figure |^ (b)) is designed so 
that the fastest growing mode is well resolved for the given resolution (256 x 256 x 512). 
The physical motivation for this initial condition is that we want to have "many" (for 
our resolution, we have around 60 small plumes at the beginning), well resolved small 
plumes seeded at the initial interface to study the mixing of the instability in details. 
Later we also show that the growth rate of this mode {k ~ 8) is close to the growth rate 
for non-diffusive, inviscid fluid (figure 

We place the position of the initial interface at z = 0.06 (—1 < z < 1), so that we can 
compare the behavior of penetrating plumes as they approach the top and bottom walls 
near the end of the simulations .|| The corresponding initial temperature profile used in 
the 2-D simulations is taken from a planar cut through the 3-D initial profile, e.g., the 
y = plane {T{x, y = 0, z)) (we also have conducted the 2-D simulations with other 
planar cuts from the 3-D profile and found no material differences). In all instances, we 
treat the case v = and n — v. 



3.1. Time evolution 

3.1.1. Mixing zone width 

We first compare the 2-D and 3-D temporal evolution; the panels in figure || show the 
2-D time evolution of the temperature, while figure ^ shows the time evolution oi y — Q 
slices from the 3-D simulations. Figure |^ provides the penetration depths as functions 
of dimensionless time for 2-D and 3-D. (To scale our time according to Youngs' time 
scale (Youngs (1991)), we multiply t with {Pg/iy^'^; thus r = {Pg/iy/'^t, where / is the 
half height of the box) . The penetration depth is determined by measuring the distance 
between boundaries of volume fraction of 1% of the hot fluid and 99% of the cold fluid. 



f From our simulation, we conclude that the effect of the top and bottom boundaries is not 
important to the propagation of the mixing zone until the mixing zone width reaches ~ 90% of 
the box size. 



384 



■11 u 



Young, Tufo, Dubey, and Rosner 
195 217 2 6 




1-= 265 300 



T= 366 370 



u u u 

r= 417 43 6 




31S 



374 




435 



3 



3 e 




4 i 



m 



1 



I 



253 



349 



•P P P P P 



396 




453 



m 



362 



405 



460 



r 



Figure 3. Time scries of the temperature field in a l{x) by 2(z) box from the 2D simulations. 
The interface between cold (black) and warm (white) fluids is perturbed by the 6T{x,y = 0) 
slice taken from the spatially random disturbances in the 3D simulations. 
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Figure 4. Time series of the temperature from the 3D simulations. The shces {T{x,v — 0, z)) 
are taken at the same times as those for the corresponding 2D frames (figure pf). 
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Figure 5. Penetration depths from 2D (a) and 3D (b) simulations for the random perturbation 
case. For the ascending structures, h is defined as the maximum value where the temperature 
reaches above 99% of the cold fluid temperature. For the descending mixing zone boundary, 
h is defined as the minimum value where the temperature drops below 99% of the hot fluid 
temperature. 

In figure ^ we show how h varies depending on different volume cuts of hot and cold 
fluids, and we conclude that penetration depths are well defined if the volume cuts 
are narrower than 1 — 99% range, which is what we have adopted to calculate results 
presented below. We plot the penetration depths against instead of r in order to 
clarify the scaling behavior. For both 2-D and 3-D, h oc ar^ during early stage of the 
evolution for 1.8 < r < 4.4. The a for this period of time is ^ 0.017 for the 2-D case, 
and a ~ 0.03 for the 3-D case. Thus, the 3-D mixing zone broadens much faster than its 
2-D counterpart. 




OuO 



^ 23 3.D ^ 4^S 




Figure 6. Penetration depths from 2D (a) and 3D (b) for different volume cuts. The 2D curves 
are from the ensemble average and show similar evolution as in figure ^(a). The diamonds are 
for 0.5 — 99.5% volume cuts of hot-cold fluids, the dashed lines are for 1 — 99% volume cuts, the 
dash-dotted lines are for 3 — 97% volume cuts, and the dash-dot-dotted lines are for 5 — 95%. 

We also observe that the scaling range 4 < < 12 is twice the range found in Youngs' 
work. The mixing zone fills the entire tank at around t = 4.7 in our case, compared 
to the values of 2.9 ^ 3.1 from simulations with low Mach numbers and high Atwood 
numbers in Youngs (1991), Finally, we observe that at time ~ 12, the increasing h 
levels off; h then resumes scaling at around ~ 16; we will discuss this phenomenon 
in detail in the following subsection. We further remark that the alpha coefficient for 
both 2 and 3-D presented above is only for early evolution. 
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3.1.2. 2-D flows versus 3-D flows: Finger geometry 

In order to understand these results, we first note that the structure of 2-D fingers 
(plumes) is horizontally stretched, whereas in 3-D, fingers remain vertical until they 
merge and form bigger plumes. This implies that the 2-D flow has a significantly larger 
horizontal velocity component when compared to 3-D; we return to this observation 
later when we discuss the evolution of single plumes. In both 2-D and 3-D, we flnd one 
large thermal pushing the envelopes at late times; however, in 3-D, we observe more 
local, small-scale structures in the interior of the mixing region than in the 2-D case. 
The large plumes push the edges of the mixing region, and dominate the expansion of 
the mixing zone at late times. The interior of the mixing zone is also very different 
between 2-D and 3-D: In 2-D, the interior structure is dominated by the large plumes, 
and therefore the general structure of the interior is simple. In 3-D, the interior is a 
combination of boundaries between plumes, residuals of mergers, and the wake behind 
the faster plumes; the interior of the 3-D mixing region is thus much more complicated 
than that of the 2-D mixing zone. A simple quantitative measure of this difference in 
structural complexity between 2 and 3-D can be constructed by measuring the stretching 
of isothermal contours during the course of R-T instability. Thus, figure ^ shows the time 
variation of the isothermal contour length (area) corresponding to T = for 2-D (3-D). 
The isothermal contour length is calculated as follows: we first identify points on the 
isothermal curves, we then calculate the length along the curve by linear interpolation 
between points. Similar technique has been utilized to calculate the iso-surface in 3D. 




Figure 7. Curvelengths of the isothermal contour (T = 0.0) for both 2D (a) and 3D (b). 



Initially, the 2-D arclength increases much faster than the corresponding 3-D arclength; 
this is because the 2-D isothermal contours are subject to horizontal stretching compara- 
ble to their vertical stretching. At late times, however, the 2-D isothermal contour breaks 
off into several connected contours, and thus a sharp decrease in the arc length results. 
In the 3-D case, the arc length increases almost monotonically until a maximum value 
is reached. This peak is due to the stretching by the shear flow as the edge plumes push 
the envelopes sideways. The break-off of this isolevel near the front explains the decrease 
in the arc length at r ~ 4.2. 

3.1.3. Energetics: 2-D versus 3-D 

We now focus on the comparison between 2-D and 3-D energetics: We flrst plot the ratio 
of kinetic energy to potential energy available to the flow as a function of time in figure 
||. Our results show a remarkable feature: the 2-D nonlinear Rayleigh- Taylor instability is 




Figure 8. Ratio of kinetic energy to gravitational potential energy: 2D versus 3D. The kinetic 
energy is the volume integral of the kinetic energy density, and the potential energy is the volume 
integral of the potential energy available in the system. For r < 1.4, the growth is exponential 
and is similar for both 2 and 3D, after this period, 2-D motions are much more efficient in 
extracting potential energy than 3-D motions. 



more efficient at extracting gravitational potential energy than its 3-D counterpart, yet 
(as shown in §3.1), the 2-D mixing zone grows more slowly than its 3-D counterpart. As we 
will show below, this app arent contradict ion is resolved when one looks at how the kinetic 
energy is partitioned. In Youngs (1991) , Youngs found the two-dimensional flow is less 
dissipative than the three-dimensional flow. In his simulations he found that the kinetic 
energy reaches a maximum at t 6 for both two and three dimensions. However, the 3- 
D flow is more "dissipative" and therefore the 3-D kinetic energy decreases significantly 
from the maximum after r = 6. He also showed that in 3-D the difference between 
potential energy and kinetic energy is greater than in 2-D. This is consistent with our 
findings. Later we show the energy spectrum of the flow inside the mixing zone from our 
3-D simulations (figure nSh. Similar energy spectrum can also be found in Youngs (1991) 




Figure 9. Early evolution of the 3D random perturbation. The diamonds are data from the 
numerical simulation. The solid line is the best fit on the log-linear plot (for 0.28 < r < 1.13), 
and the duration of the linear growth is about 4 e-folding time with a grow rate of 8. The growth 
rate obtained from the analytic results (n^ = gAk, with gA — 10 and fc ~ 8 from the power 
spectrum of our initial perturbation at the interface), is 9. 



In figure ^ we plot this energy ratio for the 3-D case at early times (the 2-D growth is 
almost the same). The growth rate from figure O is ~ 8, close to the growth rate (~ 9) for 
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fc = 8 (where the spatial perturbation spectrum peaks) from the analytic results for the 
inviscid case (Chandrasekhar (1961)). Thus one concludes that the k — 8 mode dominates 
the evolution of the initially randomly seeded perturbation. We also note that at around 
T ~ 1.3, which is 6 e-folding times after the initiation of linear growth at t '--^ 0.3, we 
begin to observe saturation due to the nonlinear terms; this marks the beginning of the 
nonlinear evolutionary stage. After this linear regime, the peak of the spectrum decreases 
in amplitude and moves toward the small wavenumber end, an indication of formation 
of large scales. 

3.2. Average Quantities and the spatial structure of the mixing zone 
3.2.1. Average Quantities 

Here we focus on the evolution of the horizontally averaged quantities: the horizon- 
tal components of the vorticities and the advective thermal fluxes uT; after horizontal 
averaging, each of these quantities is only a function of the vertical coordinate (z) and 
time (r) . We note that the buoyancy flux is the vertical component of the advective flux, 
while the buoyancy gradient is the vertical component of the advective thermal flux. 

First, we concentrate on the average horizontal components of the vorticity (figure p^ . 
The average of the horizontal vorticity is basically the vertical gradients of the horizontal 
components of the velocity, as we adopt periodic boundary conditions in the lateral 
directions. We note that the amplitude of the 2-D vorticity is almost 5 times larger than 
that of the 3-D vorticity components. The evolution of the 2-D vorticity (panel (a)) is 
indicative of the strong shear strength at the center of the mixing region. The patches 
of interchanging shades of gray in the time-space plot of the 3-D vorticity components 
indicate the merging events within the mixing region. We observe a similar pattern in 
the plot of 2-D vorticity, but only at early times; at late times, the 2-D space-time plot 
is more indicative of the roll- up of big thermals. 





Figure 10. Horizontally averaged vorticity as functions of z and t. (a) is the 2D vorticity, (b) 

is the mean horizontal 3D vorticities. 

In figure |l^ we show the time-space plot of the average advective thermal fluxes in the 
horizontal and vertical directions; the corresponding plots for 3-D are shown in figure 
|l2|. We note that in 2-D, the two fluxes are comparable in amplitude, whereas in 3- 
D, the horizontal advective fluxes are of comparable amplitude yet are both smaller 
than the vertical flux (buoyancy flux) by a factor of 5. This is a manifestation of the 
fact that in 3-D, the flow is dominantly in the preferred (i.e., vertical) direction, while 
the corresponding flow in 2-D is strongly sheared (and has comparable magnitude in 
the horizontal and vertical directions). The space-time plots for the diffusive thermal 
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fluxes basically reveal the same broadening structure of the mixing zone. The horizontal 
components in 2-D and 3-D are comparable, and the 2-D vertical diffusive flux is larger 
than the 3-D vertical flux, as found in the average advective fluxes. 




Figure 12. Horizontally averaged advective fluxes (3D) as functions of z and t. (a) is the 
mean horizontal components and (b) is the vertical (buoyancy) flux. 



To briefly summarize: We observe stronger 2-D vorticity generation compared to 3-D. 
Consequently, we observe that the 2-D horizontal transport is much more effective than 
3-D. The relative efficiency of vertical advective transport to horizontal transport, in 
contrast, is much higher in 3-D than that in 2-D. 

3.2.2. Turbulence and mixing within the Mixing Zone 

As shown in the time series for both 2 and 3-D, the spatial structure of the 2-D mixing 
zone is relatively much simpler than 3-D: the 2-D mixing zone is dominated by one or two 
big plumes while in 3-D the center of the mixing zone is full of small scale structures. This 
is clearly manifested by the temperature fluctuation probability distribution function 
(PDF) within a thin horizontal slice of the 3-D mixing zone. In our analysis we find a 
thickness of Az = 0.06 (the initial interfacial thickness) suffices to provide good statistics 
for the PDF. If we place this thin slice at the position of the initial interface, we find the 
PDF of the temperature fluctuations (AT = T — Tq) to be a gaussian centered at AT = 
(figure O). As we move this thin slice away from the initial interface, the PDF transitions 
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AT 



Figure 13. PDF of the temperature fluctuation {ST = T — Tb) within the mixing zone 
{Az = 0.025) right at the position of the original interface {z = 0.06). 



from a gaussian distribution to an exponential distribution (figure E9). Moving further 



1.0 

AT 



Figure 14. PDF of the temperature fluctuation {ST — T — To) within the mixing zone 

{Az = 0.025) a.t z = 0.21. 

than 1/3 of the mixing zone width from the center, the PDF of the thin slice transitions 
again from an exponential distribution to a yet steeper distribution (more like a delta 
function centered at either —2 or 2). The guassianity in the scalar PDF at the center of 
the mixing zone (figure |l^) implies turbulent mixing of the scalar near the center. For 
the flow inside the mixing zone, the direct evidence for turbulence is a —5/3 Kolmogorov 
energy spectrum. Due to our numerical resolution (256 by 256 by 512) the scaling range 
is small as shown in figure |l^, where the energy spectrum shows an inertial range from 
A: ~ 4 to - 18. 

As we move from the center of the mixing towards the boundary, the exponential PDF 
of the scalar at z = 0.5 implies that the flow that mixes the scalar is dominated by 
large scale. This indicates that, for the miscible RT instability and the parameters that 
we explore, turbulent mixing exists only near the initial interface, and that most of the 
mixing zone is still dominated by finger structures. These finger structures (plumes) near 
the edge are responsible for the broadening of the mixing zone. As we will show in §5, 
the propagation speed of a plume in an ambient environment is directly proportional to 
its circulation. Also shown in §5 is that, if the drag is small, these plumes initiated by the 
perturbation at the RT interface will be free falling (rising) without any turbulent scaling 
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Figure 15. Energy spectrum of the flow at the end of the simulatoin. The dashed hne is the 
Kolmogorov —5/3 scaling. The inset shows that the inertial range is from ~ 4 to fc ~ 20. 




Figure 16. Average viscous dissipation rates as functions of z and t. (a) is the dissipation 
rate due to the horizontal components of the flow, and (b) is the vertical component. 



involved. Furthermore, if the plumes pushing the edge of the mixing zone experience a 
different flow in the background, we would expect a different mixing zone broadening. 
For example, if the horizontal components of the flow are comparable to the vertical 
component, we expect the propagation of the mixing zone to slow down. Therefore, the 
transition of the ambient background flows would serve as an indication of deviation from 
the free falling phase of the mixing zone. In the next subsection we will provide evidence 
for alteration of flow patterns inside the mixing zone in the 3-D case. Combining all the 
results, we are able to demonstrate the correlation between the transition of the interior 
flows and the transition of the mixing zone broadening. 

3.3. Mixing Zone Broadening and the Internal Flow Structures 

We now focus on the relationship between the propagation of the mixing zone and the 
internal flow structures at around = 10 in the 3-D case. In figure |l^ we show the 
average horizontal and vertical components of viscous dissipation rate: the horizontal 
component is simply i/(uiV^ui -|-M2V^U2)/2 and the vertical component is (^(uaV^ua). 
In figure |l^ we show the ratio of the volume integral of these two components. 

Some observations from these results can help us gain insight into the mixing zone 
evolution: we note that average buoyancy flux (panel (b) in figure |l2|) reaches maximum 
at T ~ 2.3, while the average vertical component of the dissipation rate (panel (b) in 
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Figure 17. Ratio of the volume integrated viscous dissipation rates as a function of time (r) 
the horizontal component of the dissipation rate to the vertical component. 



figure reaches maximum around r ~ 2.8. Right after r = 2.8, we observe generation 
of vorticity inside the mixing zone, as shown in panel (b) in figure |l^. Moreover, the 
horizontal component of the advective fluxes (panel (a) in figure |2|) starts to develop 
structures around r = 2.8 as well. 

We can now summarize the evolution of the 3-D RT instability as follows: As the 
interfacial random perturbations initiates some circulation (small vortices), the seeded 
vortices increase in strength and begin free falling, accelerated by the density contrast 
across the interface (as we discuss further below). This is the free- falling phase, where the 
mixing zone is expanding due to a constant acceleration at the boundaries, and inside 
the mixing zone, the interaction between vortices is not strong enough to alter the free 
falling dynamics. At ~ 8, the dissipation rate begins to increase at a lower rate as a 
result of increasing mixing inside the mixing zone. The development of structures is most 
vigorous as the horizontal flow is enhanced and the viscous dissipation rate plateaus at 

~ 10. Hence we mark = 10 as the beginning of the mixing phase. From figure 
we also observe that the end of the mixing phase can be identified at ~ 16, when 
the viscous dissipation rate resumes positive growth. The existence of this mixing phase 
also manifests itself in the potential energy release rate and the kinetic energy partition 



(figure 18). The potential energy release rate is defined as the volume-averaged, vertical 



advective flux as follows: 



/ wTdPx = dt [ zT(fx + nT\lz\- 
Jn Jn 



(3.1) 



(We have utilized no- flux boundary conditions in deriving the above equation.) The last 
term on the right hand side is negligible for small k, as the temperature difference is 
always of order one in our simulation. In panel (a) of figure |l^, we observe just before 
the beginning of the mixing phase, the potential energy release rate reaches a plateau. 
In panel (b), the partitions of kinetic energy remain more or less constant throughout 
the mixing phase, after which they resume their original direction of varion with time. 
Thus we conclude that this mixing phase (10 < < 16) is a regime where horizontal 
components amplify and the lateral transport is also enhanced in the fiow. After this 
stage, the boundary of the mixing zone is not pushed outward by a collection of plumes 
of the smaller sizes as in the free-falling phase, but by a small number of big plumes that 
manage to survive the mixing phase and expand in size as they propagate outwards. It is 
therefore important for us to understand the evolution of plumes and their propagation 
in various kinds of ambient background flows. 
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Figure 18. Panel (a): Volume-averaged advective fluxes as a function of time. Panel (b): Par- 
tition of kinetic energy as functions of time. The solid line in panel (b) is the partition of the 
kinetic energy in the vertical direction and the dashed line is along the horizontal directions. 
The dash-dotted line is the normalized kinetic energy as a function of time. 



We also note that, based on our discussion, we expect to see departures from simple 
scaling only when the dynamics is dominated by a small number of (larger) plumes. As a 
consequence, the aspect ratio of the experimental domain can be decisive in determining 
whether such departures are actually observed: as the height-to-width ratio decreases, 
we expect to see departures from scaling at later and later times, so that plumes may 
begin to interact with the box walls before such departures can be observed. Thus we 
conclude that the dynamics and evolutions of single fingers (plumes) are important to 
the broadening of the mixing zone in both 2-D and 3-D RT instability. In the following 
sections, we first discuss RT models based on the point sources and interactions between 
each source. We will introduce the buoyancy-drag plume model in the context of the 
existent point source models, and uncover the connections between these models based 
on the physical grounds. 



4. Models for RT instability 



Early analytic work on modeling the mixing zone broadening goes back to Fermi 



(1951) 


(see also 


Sharp (1984) 


); more r 


Mukamel (1993) 


and Alon et al. (1994) 



From the numerical simulations of the random- 
perturbation cases, we observe that plumes near the edge of the mixing zone are re- 
sponsible for the broadening. We also observe that for the parameters we adopt for our 
simulations, the interior flow of the mixing zone is not homogeneous and is not isotropic. 
Instead, the flow pattern inside the mixing zone has a preferred direction associated with 
the plume structure inside the mixing zone and cannot be described as fully-developed 
turbulence. Motivated by these observations, we focus on RT models that do not require 
the existence of turbulence to achieve the free-fall scaling. We first briefly summarize two 
such models for immiscible RT instability: the first is for strongly stratified RT instabil- 
ity {A = 1)] the second is for weakly stratified RT instability (A — > 0). The similarities 
between these two models (as will be shown) lead one to conjecture that there may exist 
a point-source model that incorporates both limits. Indeed, the long-existing buoyancy- 
drag model, often used in modeling thermals or plumes, seems to capture the essential 
physical features in both models; and we will argue that this model provides an adequate 
framework for capturing the essential physics of nonlinear RT instability. We will discuss 
this third model at the end of this section. 
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4.1. 



Two models: Bubble competition model (A ~ \) 
and Point source model (A — > 0, Ag fixed). 



Zufiria's bubble competition model (Zufiria (1987)) assumes a 2-D potential flow model 
with an Atwood number of 1. In essence, this model is an extension of Wheeler's model 
( ^harp (1984) ) to the strongly stratified limit. The model consists of variables that quan- 
tify the internal potential flow within each bubble (point source inside the bubble) and 
the distances of the source of the internal flow to the bubble boundaries (radii of cur- 
vature of the bubble). For a system of 20 bubbles at the beginning, an estimated range 
0.04 < a < 0.07 is obtained from this model. 



A different model fo r 2-D inviscid, but weakly stratified RT instability, is given by Are! 
fc Tryggvason (1989) . Thi s model focuses on t he Boussinesq limit, i.e., A — > and Ag re- 
mains constant. Following Tryggvason (1988)| , the inviscid 2-D Boussinesq RT instability 
is formulated in terms of vortex sheet strength across the density interface. Vortex sheet 
strength, defined as the difference of the tangential velocity along the interface per unit 
length, is driven by the density contrast across the interface. In their model, they are able 
to capture the complicated vortex pair evolution with a simple evolutionary equation for 
the circulation per arc length. They also successfully extend the single vortex pair model 
to multi vortex pair systems, and the initial linear growth rate from their model is in fair 
agreement with the analytic results. 

4.2. Relating the models: The buoyancy-drag model 



As pointed out by Aref fc Tryggvason (1989) , their vortex model is similar to the bubble- 
competition model: the point source in the bubble model corresponds to the point vortices 
in the vortex model, and the radii of curvature of the bubble correspond to the vertical 
coordinates of the point vortices in the vortex model. The striking similarities stem from 
the fact that bubbles, especially ones elongated in the stratified direction, have flow 
patterns that amount to a vortex ring near the elongated end and a point source inside 
the bubble feeding the vortex ring. Therefore, the buoyancy-drag model used in modeling 
plumes or thermals seems to be a natural candidate to connect the bubble picture, where 
point sources are driven by a strongly stratified density contrast across the interface, to 
the vortex picture, where the point sources are driven by infinite gravitational acceleration 
(with Ag constant). 



The most simplified buoyancy-drag model (Lighthill (1986), Werne (1994)) describes 
the propagation velocity {V) of a plume as follows: 



d{ptnV) 

dt 



= Spag- CoPoutV^ 



(4.1) 



where ag is the effective acceleration due to the buoyancy contained by the thermals to 
drive the motion, Sp = pout — Pin is the density contrast between the inside and outside 
of the thermals, and Cd is the drag coefficient which depend on the geometry of the 
thermals and the fluid viscosity. There are more elaborate forms of equation (4.1), such 
as different forms for Cd, but the most essential feature is that if the correction term to 
the constant acceleration is small when compared to the acceleration, free-fall scaling is 



naturally obtained. Other variant forms of this simple model have been applied by Kane 
fe Arnett (1998) and Dimonte fc Schneider (1999)| for a variety of parameters such as 



Atwood number and Mach number, and the value of a is claimed to be satisfactorily 
in agreement with results from either experiments or numerical simulations. In the next 
section, we perturb the RT interface with small, initial vortices: Such perturbations 
generate single plumes expanding in size as they propagate away from the interface. 
As will be shown, such plumes are well described by equation (4.1). By comparing the 
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propagation velocity with the circulation of the vorticity, we are able to connect the point 
vortex model to the single plume model. We are also able to explain differences between 
2-D and 3-D miscible RT instability (either the random cases in §3 or the single plume 
cases in §5) by conducting simple energetic diagnostics on the single plumes. 



5. The evolution of plumes 

The above discussion of the buoyancy-drag model for the nonlinear development of 
the Rayleigh- Taylor instability, and the observation of Rayleigh- Taylor "fingers" in the 
experiments (§3), strongly suggests an analogy between the finger structures in the mis- 
cible version of the instability and the evolution of isolated the rmals or plume s. These 
latter structures have been investigated extensively in the past ( Turner (1973) ). Exper- 
imentally, a bubble of hot fluid is injected near the bottom of a container; the resulting 



buoyant plume expands as it rises and moves away from the source. Moses, Zocchi & 



Libchaber (1993) have explored various situations, which include injection of plumes 
into a turbulent flow. Experimental observations show that as a localized source of buoy- 
ancy is injected into the fluid, the resulting convection pattern consists of a cap which 
forms at the outwards propagating front, and a stem which connects the cap and the 
local source (Moses, Zocchi & Libchaber (1993)). In the case of the miscible RT insta- 
bility, the perturbation at the interface serves equivalently as an injection of buoyancy 
sources (or sinks) into colder (hotter) fluids. The early development of the perturbation 
then resembles a collection of plumes going up and down from the interface (figure p9[) . 



Figure 19. Plume structures at the early development of miscible Rayleigh- Taylor instability. 
Initially a layer of cold fluid (black) is on top of a layer of warm fluid (white), and a random 
perturbation is imposed at the interface. 

As shown in this figure, each rising and falling bubble has a well-deflned cap and a 
stem connecting the cap to the root at the interface, consistent with expectations for the 
behavior of plumes. 

We have numerically simulated plumes which are initiated by perturbing a Rayleigh- 
Taylor interface with a vortex pair in 2-D or a single vortex ring in 3-D. Our interest is 
primarily in plume propagation, the temporal development of the detailed geometry of 
the plumes, the dynamics and energetics of the plumes, and the dimensionality of the 
plumes. We remark in passing that 2-D plumes are referred to as planar vortex pairs in 
three dimensions, i.e., three-dimensional plumes in which all motions are restricted to lie 
in parallel, vertically-oriented planes. 

In the following, we first establish a connection between miscible RT vortex pairs and 
the immiscible version ( Tryggvason (1988)1 ) by comparing the dynamics and the detailed 
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Run number Pg v x 10^ k x 10'' d a aspect ratio 



1 


2.5 


0.5 


0.5 


0.1 


0.14 


0.2 


1 


2 


2.5 


1 


1 


0.1 


0.14 


0.2 


1 


3 


2.5 


1 


0.5 


0.1 


0.14 


0.2 


1 


4 


2.5 


0.5 


0.5 


0.1 


0.11 


0.2 


1 


5 


2.5 


1 


1 


0.1 


0.11 


0.2 


1 



Table 2. Parameters for single plume simulations; d is the interface thickness, and a is the 

amplitude of the perturbation. 



vortex structures. We then present results of two-dimensional simulations (with param- 
eters listed in Table 2), and then go on to compare these results with the corresponding 
results for three-dimensional calculations. The parameters for our single plume simu- 
lations are listed in Tabic 2. The 2-D plumes are planar, while the 3-D plumes have 
cylindrical geometry. The initial perturbation imposed at the interface is a Gaussian 
wave packet with amplitudes listed in Table 2. 

5.1. 2-D single plumes 

5.1.1. Connection to weakly stratified, immiscible R-T instability 

'iryggvason (1988) utilized the vortex sheet strength to describe the evolution of single 
mode perturbations in 2D immiscible RT instability. In the case of miscible R-T insta- 
bility, one can construct a family of isothermal contours due to the diffusivity; we can 
monitor the evolution of single plumes by evaluating various quantities (such as vortic- 
ity, thermal flux) on the isothermal (or, in Boussinesq fluids, isodensity) contours. As 



d p i 



H 1 J J 

I 

Figure 20. 2D vorticity contour plot (panel (a)) and the value of various physical quantities 
along the T = 0.75 isothermal curve (panel (b)). In panel (b), vorticity (first plot), the projection 
of the thermal gradient on the isocurve (second plot), the horizontal heat flux (third plot) and 
the vertical heat flux (fourth plot) are evaluated along the isocurve. For illustration, we place 
numbers along the isocurve to display the important spatial features related to the flow pattern. 
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long as the structure of the isothermal contours remains relatively simple (as shown in 
panel (a) of figure pO| ), we can gain some understanding of the structure of the buoyant 
plumes by considering the evolution of various physical quantities along the isothermal 
surfaces. Panel (a) shows the spatial structure of the planar vortex pair for the first 
set of parameters in Table 2. The roll-up structure is reminiscent of those observed in 



the weakly stratified, immiscible R-T instability (Aref & Tryggvason (1989)). To make 
complete comparison, we first introduce the circulation (F) along some closed curve: The 
circulation F, defined as the line integral of the flow along a closed curve, 

F= I v-dt (5.1) 
Jc 

is related to the evolution of vortex pair as follows: 

dV f f 1 , , 

where A is the surface enclosed by the curve C . In the case of weakly stratified Rayleigh- 



Taylor instability ( Aref fc Tryggvason (1989) ), the above relation is modified by replacing 



the right-hand-side integral with the potential energy density due to the interfacial den- 
sity contrast. The propagation velocity {V) of the vortex pair in a nonstratified fluid is 
then related to the circulation F as 

V = F/77, (5.3) 

where 77 is basically the correlation length for vortices, and is a function of geometry 
and dimensionality: one expects 77 for 2-D vortex pairs to be larger than for 3-D vortex 



rings. Equation (5.3) is valid even in the case of weakly stratified R-T instability, where 
vortices, originally initiated at the interface due to the density perturbation, propagate 
in a uniform ambient background. 

With equation ( |5.3| ), we can connect the bubble model to the point vortex model via 
the buoyancy-drag plume model. To fix ideas, we first outline an evolutionary scenario 
for the weakly stratified, miscible R-T instability: As we disturb the interface with some 
temperature perturbation, vortices are generated as a result of lateral thermal gradient. 
Thus some initial circulation is seeded inside the vortices and a propagation velocity is 
acquired according to equation (^.3|) . The circulation then evolves according to equation 



(5.2), and the propagation velocity is approximately given by equation (5.3) at each 



time step for an updated circulation (Hill (1975)). We also observe that the acceleration 



of circulation F. Combining equations ( ^.2| ) and (5.2), we obtain the equation for the 
acceleration for the planar vortex pair: 



experienced by the vortex pair is proportional to the increase/decrease rate of change 

; equations ( ^ ) 
mrtex pair: 

dV d(T/r]) 



dt dt 



gfiSp,s), (5.4) 



where g in equation (^) is the gravitational acceleration, and / is a function of density 
contrast {Sp) and geometry of the interface (unit tangent vector along the interface s). 
For point vortices ( Aref fc Tryggvason (1989)| ), f = g ■ s where g is the unit vector 



along the gravitational direciton and s is the unit tangent vector along the interface. In 
our case, we expect / to be a constant (effective acceleration) plus some correction term 
due to viscous drag. This expectation rises from two facts: The first is the observation 
that the plume consists of a vortex pair inside the cap and a flow sustaining the vortex 
structure (the neck connecting the cap to the stem). Therefore, we expect the dynamics 
of the vortex pair to be similar to the plume dynamics. The second fact is that the 
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Figure 21. Circulation T (along the T = 0.7 isothermal contour) and the propagation of a 
2-D single plume. In panel (a), the solid line is from the simulation, and the dashed line is the 
fit. In panel (b), the solid line is from the simulation, the dashed line is the fit, and the 
dash-dotted line is the fit to the single plume model. 



propagation of buoyant plumes (not sustained by a constant source) in a homogeneous 



fluid is described by the buoyancy-drag model (equation 4.1), where the constant density 
difference across the plume boundary provides buoyancy acceleration, and the viscous 
drag experienced by the plume causes a deceleration. To be more specific, we expect / 
to take the following form: 

f = a + ji,y,Sp,V), (5.5) 

where the "constant" a may depend on the dimensionality and the initial perturbation, 
and the correction terms are included in the function 7 (a function of viscosity ly, the 
propagation velocity V and the density difference dp) . If 7 is chosen to be proportional to 
—V^ as in Werne (1994)| , the propagation velocity will reach a constant terminal velocity 
{V — Vf is a. constant) after a free falling period {V oc t and h (xt^). 

To fully establish the correspondence between the immiscible point vortex model and 
the miscible plume model, we calculate the circulation F along isothermal contours (not 
necessarily closed) in our 2-D simulations. Figure (a) shows the temporal evolution 
of F along T = 0.7, and during 0.7 < r < 2.1, we find F to be proportional to r^. 
Comparing with figure pll (b) , we note that this duration coincides with the free falling 



phase when h scales to r . This is similar to what is observed in Aref & Tryggvason 



(1989) : the propagation of the 2-D vortex pair is free-falling at the beginning, and it 
decelerates due to the decrease of available potential energy per unit arclength as the 
rollup continues (and also due to diffusion and viscous dissipation in our case) - thus 
the decrease in circulation inside the vortex pair. We also note that, the constant initial 
acceleration rate of F even in the miscible case indicates that the initial circulation is 
due to the buoyancy-driven flow around the cap head. As the circulation increases, drag 
becomes important and the circulation saturates to a maximum. The circulation then 
decreases at late times due to the strong enhancement of the horizontal motion, which is 
typical of 2-D flow and less so in 3-D flow (we will discuss this difference later) . We should 
point out here that, despite the similarities already presented between point- vortex model 
and plume model, there exists an essential difference: the plume has a stem structure 
associated with the circulation around the neck, which closely resembles the flow in the 
vicinity of RT "flngers" . Our future research problem is to understand how the presence 
of multiple plumes affects the circulation near the necks; this will be an important future 
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step in connecting these models to the RT finger evolution in more vigorous flows (higher 
Reynolds numbers). 

We are now ready to compare the vortex structure between the miscible R-T insta- 
bility and the immiscible version. From panel (b) of figure ^ we obtain the following 
observations: 

(1) The vorticity vanishes at the tip of the cap head. This implies that the flow near the 
cap is essentially irrotational, and explains why the potential flow assumption adopted in 



Moses, Zocchi & Libchaber 



the der ivation of the cap shape by Lighthill (1986) (see also 
(1993) ) works so well in describing the head sha pe. This is also consistent with results 
for the weakly stratified, inviscid R-T instability ( [Aref fc Tryggvason (1989)1 , Tryggvason 
(1988)1) 



(2) The vorticity reaches extrema near the neck of the plumes, i.e., near the part of the 
plume which connects the stem and the cap head; this is the location where circulation 
is most vigorous. This is also consistent with the fact that circulation reaches extremum 
at the inflection points along the interface ( Aref fc Tryggvason (1989) ). 

(3) The amplitude of the thermal gradient reaches extrema either near the cap heads or 
near the neck of the plumes. This is consistent with the previous observation, because the 
larger the thermal gradient is, the larger the rate of circulation is, and thus the vorticity 
reaches extrema at these points. 

5.1.2. Dynamics and energetics of 2-D plumes 

In this subsection, we study the details of the single plume structures to uncover the 
physical ingredients essential to model the RT fingers. As in the point vortex model, we 
find the evolution of arclength of isothermal contour to be indicative of the detailed evolu- 
tion of the plume structures. Here we first focus on the evolution of the isolevels associated 
with single plumes. Secondly we explore the dependence of the single plume evolution on 
the various parameters such as the amplitude of the perturbation, the Reynolds number 
and the Prandtl number. 
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Figure 22. Amplification of arc length as a function of time for 15 different isolevels (2D). 

Figure ^ is an example of how the arc length evolves as a function of time for different 
levels from the 2-D simulation (run number 4 in Table 2). These are the time evolution of 
arc length for 15 different isothermal contours. These various evolutions in the isothermal 
contours can shed light on the detailed, local structure of the fiuid. First we observe that, 
among all the 15 curves, curve 11 (T = —0.4) has the longest curve length. Curves 1 — 6 in 
figure ^ correspond to isolevels near the boundary between the plume and the ambient 
fiuids. Before t ~ 1.0, the fiow inside the plume pushes these levels upward as the 
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plume rises. From r ~ 1.0 on, these contour levels begin to catch up with the ascending 
boundary of the plume to the heavy fluid, and the level sets break up into several smaller 
connected loops. After this event (at r ^ 1.4), the plume has established its own shape 
and the cap propagates in a shape-preserving fashion. At time 1.7 < r < 2.1, small 
vortex pairs near the original interface begin to roll up to bigger vortex pairs and thus 
the increase in arc length. Curves 7 — 11 represent transition contours inside the plume. 
The sudden drop at r ~ 2.1 correspond to the shedding of smaller vortex pairs. Curves 
12 — 15 represent the evolution of the levels near the boundary of the ascending plumes, 
and peaks in each curve represent break-offs of isolevel curves. 

We now focus on how the penetration depth h grows in time, and how this temporal 
evolution depends on various control parameters in the simulations. Figure ^ shows 
curves of h versus r for our 2-D simulations. We first note that at early time, the ascending 
h (top panel of figure [2^) does scale with from ~ 0.2 to ~ 3.0. After = 3.0, the 
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Figure 23. Penetration depth versus for the 2D single plumes. The top panel is for the 
ascending boundary and the bottom panel is for the descending boundary. The numbers are the 
run numbers listed in Table 2. 

plumes have reached a terminal velocity and therefore travel at a constant speed, which 
is often observed in plumes without sustaining sources. We further note that the higher 
the Reynolds number, the faster the plume penetrates (cf. curve 1 and curve 2). We 
also observe faster penetration as we increase the diffusivity ratio (y I ti) (cf. curve 2 and 
curve 3), which is a direct result of weaker stabilization for smaller diffusivity. We have 
also varied the initial perturbation amplitude (cf. curves 1 and 4, and curves 2 and 5): 
For perturbation of amplitude below 8% of the total box height, we find that the larger 
the initial perturbation, the faster the plume penetrate. The penetration depth depends 
linearly on the initial amplitude until r ~ 1.4, where nonlinearity may have set in to 
amplify the dependence on the initial state (cf. curve 2 and curve 5). We also note that the 
penetration of the heavy fluid into the light fluid remains more or less similar as we vary 
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these parameters. The strong dependence of the ascending structure on the parameters is 
not reflected on the dynamics of the descending structures. This demonstrates that, with 
this form of initial perturbation, the ascending structures have different dynamics than 
the descending structures despite the low Atwood number adopted under the Boussinesq 
approximation. This is certainly not the case if we perturb the interface with symmetric 
disturbances. 



5.2. Comparison between 2-D and 3-D plumes 




Figure 24. Comparison of 2D (panel (a)) and 3D plume (panel (b)) structures: geometric 
shape. The solid line is the T = 0.7 isothermal contour. 

In this section we present results from simulations of both 2-D and 3-D single plumes 
for the first set of parameters in Table 2. The differences between 2-D and 3-D plumes 
are manifested in various aspects. In figure ^ we display the plume structures (at the 
same penetration depth) in both 2-D and 3-D. We also compute various quantities along 
the isothermal contours, and basically the same spatial structures are found in the 3- 
D plumes: the vorticity goes to zero at the cap head, the maximum circulation occurs 
around the neck, and positions of where temperature gradient have extrema are also 
points where the vorticity have extrema. The only discernible difference between 2-D and 
3-D is the magnitude: the 3-D thermal gradient has higher maximum at the cap head, 
and the advective heat transport (both vertical and horizontal components) in 3-D is 
much smaller compared to the 2-D heat transport. The latter (difference in the efhciency 
of heat transport) explains the thinning of the 2-D stem as opposed to a constant width 
of the 3-D stem observed in our simulations. From our simulations, we also conclude 
that, for the same parameters, the cylindrically symmetric plumes (our 3-D plumes) 
travel much faster than the planner plumes (our 2-D plumes). As seen in figure |2^, the 
penetration depth in 3-D increases faster with time than in the 2-D case, and eventually 
reaches almost 2 times the 2-D penetration depth at the end our our simulations. In order 
to understand this difference, we have calculated the temporal evolution of the ratio of 
total kinetic energy to the potential energy, and the partition of kinetic energy (which 
we break up into its two components, namely the total kinetic energy associated with 




Figure 26. Panel (a): Ratio of total kinetic energy to total potential energy available in the 
system as a function of r. Panel (b): Vertical partition of the total kinetic energy and horizontal 
partition of the kinetic energy as functions of r. In both panels, we plot 2D (dashed lines) 
together with 3D (solid lines). 



vertical motions and with horizontal motions). Thus, in figure (a) we plot the energetic 
ratio versus t, and show the evolution of the kinetic energy partition as a function of 
time in figure (b). First, we point out that the rate of conversion of potential energy 
to kinetic energy is roughly the same in 2 and 3-D. Secondly, we note that partitioning of 
the kinetic energy into horizontal and vertical components is different in 2 and 3-D (with 
a larger fraction of the kinetic energy going to horizontal motions in the 2-D case). Based 
on these facts, we can now attribute the increased effectiveness of plume penetration in 
three dimensions to the fact that while more potential energy is released per unit time 
in 2-D than in 3-D, the partitioning of the resulting kinetic energy into its components 
drastically favors the directed vertical component in 3-D, but not in 2-D. into the vertical 
component of the flow kinetic energy, so that the 3-D plumes go faster and penetrate 
deeper in the same amount of time than their 2-D counterparts. 

5.3. Application of the single-plume model and future direction 

We apply the plume model to our single plume simulations and find "good" agreement 
as we fit the penetration depth h to the model (in the sense that the absolute error is less 
than 1% over the whole evolution). A typical example is shown in figure |T](b), where the 
dashed line is the least square fit for the free-falling law, and the dash-dotted line is the 
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least square fit for the plume model. We clearly see that the plume model captures the 
whole evolution, and the deviation is small compared to the fit. When we apply the 




Figure 27. The 3D random case revisited: the diamonds are from the simulations and the 
dashed lines are the fits of the analytic plume model to the data during the first free-falling 
phase. 



^ < 13.4, where the effective acceleration 
n, the evolution deviates from the fit at 



plume model to the random case 27, the existence of three evolutionary phases prohibits 
a physically meaningful fit for the whole duration of the evolution, and a reasonably 
good fit is found for early evolution: 2.6 < r 
from this fit is a ^ 0.03. As shown in figure 

~ 11, right after the beginning of the mixing phase as discussed in §3. This is an 
indication that the interaction of plumes with the background fiow is substantial enough 
to arrest the free falling propagation of the mixing zone. To incorporate this interaction 
with the background fiow into the single plume model, one needs an equation for the 
energy evolution of the plume, and also an additional momentum term in equation 4.1 to 



describe the momentum transfer between the horizontal motion and the vertical motion. 
The interaction between plumes, such as merging, also has important contribution to the 
evolution of plumes as they expand in size and getting closer to each other. This is a 
challenging future direction and is currently under investigation. 



6. Summary and Conclusions 

We have performed numerical simulations in both 2-D and 3-D with relevant length 
scales arranged such that the thickness of the interface between the two fiuids remains 
the same throughout the simulations. We have found dramatic differences between 2-D 
and 3-D flows with the aid of average quantities. 

We first perturb the RT interface with "random" perturbations and investigate the 
evolution of the mixing zone in both 2 and 3-D with the aid of averaged quantities. 
Scaling of the mixing zone width with is found in both 2 and 3-D, but only during 
selected intervals of time. Moreover, we find the 3-D mixing zone expands two times 
faster than the 2-D mixing zone. An important finding in this respect is that we can 
identify three phases of evolution for the 3-D mixing zone: the first is the free-falling 
phase right after the linear growth period (after 6 e- folding times, §3), the second is 
the mixing phase where the free falling slows down and more mixing is generated as a 
result of enhanced horizontal motions, as illustrated in §3; the third is another free falling 
phase during which the broading of the mixing zone resumes the scaling, when the big 
plumes near the edge become decoupled from the mixing zone and propagate on their 
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own. Whether aU of these phases are observed is hkely to depend on the aspect ratio 
of the experimental domain. In the 2-D single plume case, we explore the parameters 
to demonstrate the effect of initial conditions on the broadening of the mixing zone and 
thus the effective acceleration a. Through energetic analysis, we are able to attribute 
the difference between 2 and 3-D flows to the fact that in 2-D, more energy goes into 
highly correlated horizontal and vertical motion (i.e., vortical motions) than into directed 
motions along the gravitational direction (3-D). 

We have discussed two analytic Rayleigh- Taylor models relevant to our computations 



(the bubble-competition model (Zufiria (1987)) and the point-vortex model (Aref & 
Tryggvason (1989)] ), which arc effectively two extreme limits of the 2-D incompressible 
RT instability from the point-source point of view. The essential physical features found 
in both models are also key constituents of the buoyant-drag plume model, in which 
each single plume experiences an acceleration due to the density contrast, and also a 
deceleration, which is due to a combination of viscous drag and the horizontal shear 
flows in the ambient background. We discuss how the plume model serves to connect the 
two models on physical grounds, and point out that the plume model actually contains 
more physical features, such as the coupling of the plumes with the mixing zone via the 
interaction between a background flow (inside the mixing zone) and flows in the stems of 
the plumes (whose heads are pushing the envelope of the mixing zone). Furthermore, we 
show that for the miscible vortex pair, both circulation (F) and penetration depth (ft.) 
scale with for the same duration of time. This illustrates the close connection between 
the plumes and the point vortex pairs. In addition, we also strengthen the connection by 
comparing the evolution of RT interface in the miscible case with that in the immiscible 
case. The extension of the single plume model to multi-plume systems is of great interest 
and is now under investigation. 

The success of the buoyant-drag model is an indication of the fundamental physical 
ingredients of the RT instability. Each bubble can be considered as an internal point 
source, and the curvature of the bubble evolves according to the internal source strength 
and the relative position of the internal source to the bubble boundary. The internal point 
source is fed by the flow from below the bubble, namely the stem connecting the cap 
head of the bubble. The coupling of the stem with the background flow, then, determines 
the strength of the internal source. All these features are ingredients of the plume model, 
which is, strictly speaking, an empirical model for plumes in both laminar and turbulent 
flows. Therefore, we find the buoyant-drag plume model to give an adequate description of 
the various physical processes that we have observed. The extension of this single plume 
model to incorporate the coupling of plumes, merger events, and the background flow 
interaction is challenging and will no doubt bring us more insight into the RT instability. 
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